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Abstract 

Bose-condensed systems with broken global gauge symmetry are considered. The 
description of these systems, as has been shown by Hohenberg and Martin, possesses an 
internal inconsistency, resulting in either nonconserving theories or yielding an unphys- 
ical gap in the spectrum. The general notion of representative statistical ensembles is 
formulated for arbitrary statistical systems, equilibrium or not. The principal idea of 
this notion is the necessity of taking into account all imposed conditions that uniquely 
define the given statistical system. Employing such a representative ensemble for Bose- 
condensed systems removes all paradoxes, yielding a completely self-consistent theory, 
both conserving and gapless in any approximation. This is illustrated for an equilib- 
rium uniform Bose system. 
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1 Introduction and Analysis of Problem 



The appearance of Bose-Einstein condensate in a Bose system is usually associated with 
the breaking of the global gauge symmetry, which is commonly realized by means of the 
Bogolubov shift in the field operators [1-4] . The idea of symmetry breaking is in line with 
the general understanding that phase transitions of different nature are accompanied by some 
symmetry changes. In theoretical description, there exist several ways of symmetry breaking 
(see review article [5]). The most known is the Bogolubov method of quasiaverages, when 
the symmetry of the Hamiltonian is disturbed by infinitesimal sources [4] . For breaking the 
global gauge symmetry, the latter technique is not always convenient, while the Bogolubov 
shift is a simple and sufficient condition for the symmetry breaking [6]. 

When the global gauge symmetry is broken, then the description of Bose systems, based 
on the standard grand canonical ensemble, encounters a dilemma, first discussed by Hohen- 
berg and Martin [7], who emphasized that this description results in either nonconserving 
theories or yields an unphysical gap in the spectrum of particles. In nonconserving theories, 
local conservation laws are not valid, which, at the same time, is connected with inconsis- 
tencies in thermodynamics. 

The Hohenberg-Martin dilemma of conserving versus gapless theories is usually illus- 
trated by considering some approximate schemes [7], for instance, a kind of mean- field ap- 
proximations. The origin of this dilemma can be explained as follows. 

Bose-Einstein condensation implies the macroscopic occupation of the ground-state en- 
ergy level, when the number of condensed particles Nq is such that the condensate fraction 
Nq/N is nonzero in the thermodynamic limit. The total number of particles in the system, 
N — Nq + Ni, becomes a sum of the condensed-particle number A^o and the number A^i of 
uncondensed particles. According to Bogolubov [1-4], the number of condensed particles A^o 
is such that it provides thermodynamic stability for the system, minimizing the thermody- 
namic potential. At the same time, Nq = N — Ni, where A^i = Ni{fi,T, p) is a function 
of the chemical potential /x, temperature T, and density p = N/V, with V being the sys- 
tem volume. Hence, A^o = Nq{p,T,p) is a function of the same variables. Conversely, the 
chemical potential // = //(T, p) is a function of T and p. Thus, the condition of the stability, 
realized as the minimization of the thermodynamic potential, defines the chemical potential 
p. It is important to emphasize that this procedure of minimization does not depend on 
approximations, but is generally valid, as was strictly proved by Ginibre [8]. 

Prom another side, the particle spectrum, under the broken gauge symmetry, has to be 
gapless. This is, actually, the condition for the existence of stable Bose-Einstein condensate. 
Since, if there would be a gap in the spectrum, there could be no macroscopic occupation 
of a single ground-state level. Hugenholtz and Pines [9] showed that the chemical potential 
p = Ell (0,0) — 5]i2(0,0) is expressed through the self-energies EQ^(k, a;), and this relation 
makes the spectrum gapless for any Bose system. As far as the self-energies Sq,^(0,0) are 
the functions of T and p, the Hugenholtz- Pines relation defines the chemical potential p — 
p{T, p). 

In this way, the condition of thermodynamic stability, in the frame of the Bogolubov- 
Ginibre minimization procedure, defines a chemical potential that we shall denote as Pbg- 
At the same time, the existence of a gapless spectrum, connected to the Hugenholtz- Pines 
relation, also defines a chemical potential, which may be denoted as php- These chemical 
potentials, found in two different ways do not necessarily coincide, and, generally, they 
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are different: ^bg 7^ A*i?p- Therefore, if one accepts as tlie system diemical potential 
tlie Bogolubov-Ginibre value fiBG^ providing thermodynamic stability, then one acquires an 
unphysical gap in the spectrum, proportional to the difference \iihp ^ f-i'Scl^^'^- Conversely, 
accepting as the chemical potential the Hugenholtz-Pines form /j-hp, one gets a nonconserving 
theory with incorrect thermodynamics, since the stability condition does not hold. This is 
the origin of the Hohenberg-Martin dilemma of conserving versus gapless theories [7]. 

In any case, the fact that fJ^sG f^HP makes the system unstable and its description 
not self-consistent. This problem docs not arise only in the limit of asymptotically weak 
interactions, when the Bogolubov weakly- noni deal gas approximation is applicable [1,2]. In 
this limit, /ibg and /ihp asymptotically coincide. However, in any higher approximation one 
has II BG (J'HP- 

For practical applications, one usually does the following. When one is interested solely in 
the system dynamics, but not in its spectrum, one derives the equations from a variational 
principle, which guarantees the validity of conservation laws [10,11]. For an equilibrium 
system, this is equivalent to the Bogolubov-Ginibre variational procedure. And, when one 
studies only the system spectrum, one accepts the Hugenholtz-Pines relation, forgetting 
about inconsistent thermodynamics and instability. Clearly, such palliative ways are not 
satisfactory. The principal problem remains how to make the general theory both conserving 
as well as gapless. 

There have been several attempts to cure the problem, which could be classified into 
three groups: 

The most often used trick is the omission of anomalous averages. As is clear, this is a 
rather unjustified way, since, as soon as the global gauge symmetry is broken, the anomalous 
averages do exist and are not zero. One, sometimes, ascribes this trick to Popov, calling 
it the "Popov approximation". However, it is sufficient to look attentively at the original 
works by Popov [12-14], which are usually cited in this respect, in order to realize that he 
has never suggested anything like that. He considered the properties of a Bose gas in the 
vicinity of the critical temperature T^, honestly calculating all terms, normal and anoma- 
lous. When temperature tends to T^,, the condensate density tends to zero together with 
the anomalous averages. Actually, both these quantities, the condensate density and the 
anomalous averages are the order parameters, appearing together in the broken-symmetry 
phase, and disappearing also together, when the symmetry is getting restored at the critical 
temperature. Contrary to this, the normal average, that is, the density of uncondensed par- 
ticles, increases when approaching the critical temperature, reaching at Tc the total system 
density. This is why at the close vicinity of T^, the anomalous averages become, without any 
special assumptions, smaller than the normal ones. However at low temperatures T <S Tc, 
the anomalous averages not merely can be of the order of the normal averages, but can even 
be much larger than the latter, as direct calculations show [15]. Moreover, omitting the 
anomalous averages makes the system principally unstable [15,16]. Thus, this trick neither 
has anything to do with Popov nor can be accepted as a reasonable approximation at low 
temperatures. Additionally, if one would wish to ascribe a name to this trick, one should 
know that the first person, who really suggested and used it, was Shohno [17]. And this 
was known in literature as the Shohno model. The word "model" is appropriate here, since 
this is, actually, just a model, but not a justified approximation. For example, Reatto and 
Straley [18] used the term "Shohno model" and studied its properties. 

Another way of removing the gap in the spectrum is to calculate the chemical potential 
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and self energy in different approximations. Then one defines the chemical potential from 
the thermodynamic stability condition in one approximation, but for calculating the self- 
energy, one invokes a higher-order approximation, so that to cancel the gap. The additional 
higher-order terms can be motivated by Bethe-Salpctcr or scattering-matrix approximations 
[19,20]. This way is what was called by Bogolubov [4] the mismatch of approximations. 
Bogolubov already mentioned [4] that such a mismatch can really influence the appearance 
or disappearance of the spectrum gap, but renders the theory not self-consistent and the 
system unstable. 

One can also kill the spectrum gap by adding to the self-energy phenomenological terms, 
such that to cancel the gap [21-23]. This way is, clearly, equivalent to the previous one, 
since changing the self-energy by invoking some higher-order approximations is not unique 
and, hence, is also phenomenological. 

As is evident, all these ways, attempting to cure the problem, are based on a kind of the 
mismatch of approximations and, as a result, they have the same common defects: 

(i) They are not uniquely defined since there exists an infinite number of particular tricks 
for removing the spectrum gap. So, the ambiguity remains [23]. 

(ii) They are not self-consistent, involving in this or that way the mismatch of approxi- 
mations [4]. 

(iii) They, as a rule, correspond to an unstable system, either with a not minimal ther- 
modynamic potential or with a divergent susceptibility [15,16]. 

(iv) The order of the condensation transitions changes, resulting in a first-order phase 
transition, instead of the correct second-order one. This happens because of the internal 
inconsistencies in the description. The disruption of the phase-transition order is a common 
feature of inconsistent approximations, which either do not satisfy the stability conditions 
or possess a spectrum gap [18,21-23]. This was already noticed in the early works analysing 
the Hartree-Fock-Bogolubov approximation having a gap in the spectrum [18,24-27]. The 
generality of such a change of the phase-transition order from second to first in different 
mean-field approximations was emphasized in a detailed discussion by Baym and Grinstein 
[28] and recently by Kita [23]. As is clear, the thermodynamics of a system with a first- 
order phase transition is rather different from that of a system displaying a second-order 
transition. And in the vicinity of the critical point, the thermodynamics in these two cases 
differs drastically [23,28]. 

(v) Finally, the mismatch of approximations, by its own, is not a regular procedure. For 
each given approximation, it is necessary to invent special tricks, which, as is stressed above, 
are not uniquely defined, hence, ambiguous. There is no a general rule how to do this in a 
unique way for approximations of different order. 

Thus, we have to conclude that the Hohenberg-Martin dilemma [7] remains unsolved. 
The methods, based on the mismatch of approximations, are not self-consistent, having 
several internal defects discussed above. In order that an effective theory be self-consistent, 
all dynamic and thermodynamic equations must be derived from the same Hamiltonian, or 
Lagrangian, and treated in one chosen approximation, without involving the approximation 
mismatch [29,30]. The most recent thorough discussion of the Hohenberg-Martin dilemma, 
with many citations, can be found in the review article by Andersen [31]. 

In the paper [6], the idea was advanced that the problem can be solved by employing a 
representative statistical ensemble. In Refs. [32,33], it was shown how to make the Hartree- 
Fock-Bogolubov (HFB) approximation for a dilute gas both conserving and gapless. The 
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aim of the present paper is to develop a general approach, independent of particular ap- 
proximations, for the self-consistent treatment of arbitrary Bose systems with broken gauge 
symmetry and to demonstrate on the most general footing that the resulting theory is really 
completely self-consistent. The approach is based on the notion of representative statisti- 
cal ensembles, whose general formulation is given in Sec. 2 for both equilibrium as well as 
nonequilibrium systems. This notion is specified in Sec. 3 for Bose systems with broken 
gauge symmetry. Thermodynamic self-consistency of the theory is emphasized in Sec. 4. 
Operator equations of motion arc obtained in Sec. 5 for a Hamiltonian with an arbitrary 
interaction potential. Sec. 6 demonstrates that the local conservation laws are valid on the 
operator level, hence, being automatically satisfied for the related average quantities. The 
equation for the condensate wave function is derived and analyzed in Sec. 7 for an arbitrary 
Bose system. In Sec. 8 a uniform Bose system is considered and illustrated for the HFB 
approximation. The behaviour of the condensate and superfluid fractions is studied in Sec. 
9. The equations for Green functions are presented in Sec. 10, where it is shown, by using 
the Bogolubov theorem [4], that the Hugenholtz-Pines relation follows, thus, proving the 
complete self-consistency of the developed approach. Sec. 11 is the conclusion. 

Throughout the paper, the system of units is used, where the Planck and Boltzmann 
constants are set to unity, h = 1, ks — 1- 

2 Representative Statistical Ensembles 

The idea of representative statistical ensembles goes back to Gibbs [34], who mentioned 
that to prescribe a canonical distribution for a system may be not sufficient, but this distri- 
bution has to be complimented by those constraints and conditions that provide a correct 
representation of the considered statistical system. The term "representative ensembles" 
was used by ter Haar [35,36], who investigated the problem of a proper representation of 
equilibrium statistical systems. Equilibrium and quasiequilibrium representative ensembles 
were described in the review article [5]. In this section, we shall formulate the notion of 
representative statistical ensembles for both equilibrium and nonequilibrium systems. 

Representative statistical ensembles for equilibrium systems are sometimes also termed 
as generalized Gibbs ensembles, subjective or conditional ensembles. Their mathematical 
construction is based on the conditional maximization of the Gibbs entropy, as was done by 
Janes [37,38]. This concept was also employed by Girardeau [39]. 

Let us start with equilibrium systems. An equilibrium, stafAstical ensemble, by definition, 
is a pair {JF, p} of the space JF of microstates and a statistical operator p. In order to 
correctly define the latter, it is necessary, according to Gibbs [34], to take into account all 
conditions and constraints, imposed on the system. Suppose, we have a set {Q} of self- 
adjoint operators, defined on the space which will be called condition operators. This is 
because these operators enter the statistical conditions 

C, = < (7, > = TVp Q , (1) 

which are necessary to take into account for correctly representing the considered system. 
The trace operation in Eq. (1) is over the given space !F, which can be defined as a Fock 
space. The first evident constraint is the normalization condition for the statistical operator, 

1 = < > = Trp , (2) 
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where Ijr is the unity operator in The definition of the internal energy 



E = < H > = Tip H , (3) 

as the average of a Hamiltonian H, is another standard statistical condition. But, in addition 
to constraints (2) and (3), there can exist any number of other statistical conditions (1). The 
conditional maximization of the Gibbs entropy 

S'=-Trplnp (4) 

is equivalent to the unconditional minimization of the information functional [40] , defined 

I[p]^-S + Xo{Trp-l)+(3(TrpH-E)+(3Y.''i{^pCi-Q) , (5) 



as 



in which the standard conditions (2) and (3) are included exphcitly. The quantities Aq, P, 
and are the appropriate Lagrange multipliers, with /3 = 1/T being inverse temperature. 
Minimizing the information functional (5) with respect to the statistical operator p gives 

P - \ e-^- , (6) 
where Z = exp(l + Aq) is the partition function and the grand Hamiltonian 

H = H + Y.^A (7) 

i 

is introduced. The representative statistical ensemble, under constraints (1), (2), and (3), is 
then the pair {T^p} of a space T of microstates and the statistical operator (6), with the 
grand Hamiltonian (7). When one of the condition operators Ci is the number-of-particle 
operator N and the related Lagrange multiplier i/j = — one gets a particular form of the 
grand Hamiltonian H = H — p,N. However, any other necessary constraints can be included, 
resulting in the general expression (7) for the grand Hamiltonian. 

The condition operators Ci are to be self-adjoint, Cf = Ci, so that the grand Hamiltonian 
(7) be also self-adjoint. In many cases, the condition operators are taken as integrals of 
motion, such that [Ci, H] — 0. But this is not compulsory. For instance, the number-of- 
particle operator N does not commute with the Hamiltonian energy operator H, when the 
global gauge symmetry is broken. 

The representative ensemble {J-',p}, with the statistical operator (6) and the grand 
Hamiltonian (7), define all thermodynamics of an equilibrium, or stationary, system. The 
construction of this ensemble has been more or less straightforward, following the ideas of 
Gibbs [34], ter Haar [35,36], and Janes [37,38], as is reviewed in Refs. [5,40]. But the defini- 
tion of representative ensembles for arbitrary nonequilibrium systems is not evident. Below, 
we give the generalization of the notion of representative ensembles for nonequilibrium sta- 
tistical systems. 

To describe a nonequilibrium system, one needs, in addition to the space of microstates 
JF and the initial value of the statistical operator p = p(0), to define the temporal evolution 
of the system. This evolution can be described by the time-dependent statistical operators 
p{t), satisfying the Liouville equation, or by the time dependence of physical operators, sat- 
isfying the Heisenberg equation. Equivalently, the time evolution can be associated with the 
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evolution operator, satisfying the Schrodinger equation. Keeping in mind any of these ways, 
we may denote the prescribed temporal evolution by the symbol d/dt. Then a nonequilih- 
rium statistical ensemble is a triplet {JF, p, d/dt}. Clearly, when the time evolution is absent, 
or trivial, this definition reduces to that for the equilibrium case. 

To be more specific, we may remember that each system is characterized by some dy- 
namical variables, such as field operators. Let us keep in mind a set of field operators ip{x, t), 
whose particular representation is not important at this stage. For example, the variable x 
can represent real-space coordinates or momenta. It may also include other continuous or 
discrete variables, such as the spin indices or component labels. All physical operators, such 
as the Hamiltonian energy operator Hlip], are functionals of the field operators. 

The most general way for describing the system dynamics, as is known [41], is the ex- 
tremization of the action functional. Implementing this for our case, we need the Lagrangian 



in which the notation 



L[i;] = E[ij]-H[^p], (8) 

Eit/j] = J iP\x, t) i ^ iP{x, t) dx (9) 

for the temporal energy operator is used. 

To make the ensemble representative, we have to take account of all additional conditions 
uniquely characterizing the system. This implies that, similarly to Eq. (1), we have to take 
care of the statistical conditions 

C, =< Ci[i,] > , (10) 
where C'Jt/'] are the appropriate condition operators and 

< Ci[iP] > = TVp(O) Ci[iP{x,t)] = Trp{t) Ci[iP{x,0)] . 

The principle of action extremization, under the given statistical conditions (10), is equiv- 
alent to the unconditional extremization of the effective action 

Aiij]^ I {Liij] - i^Am} dt (11) 

with the Lagrange multipliers Ui guaranteeing the validity of conditions (10). Combining 
Eqs. (8) and (11), we can rewrite the effective action (11) as 

A[^]^ J {m-H[^]} dt, (12) 

with the grand Hamiltonian 

H[^l;]=H[^l;] + Y:i^iCi[^], (13) 

i 

having the same form as in Eq. (7). The extremization of the action functional with respect 
to field operators means the variational equation 

= (14) 
5jlj^{x,t) ^ ' 
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plus its Hermitian conjugate. Equation (14), in view of action (12), is identical to the 
equation 

■ d , SHli/j] 

The evolution equation (15) is what one needs for a complete definition of the nonequilib- 
rium representative ensemble. It is important to stress that the system dynamics is governed 
by the same grand Hamiltonian as its thermodynamics. 

We may also note that in the Heisenberg representation, as is well known (see, e.g., Refs. 
[40,41]), the variational equation (15) is the same as the Heisenberg equation 

i^^^{x,t) = [^(a;,t), H[ilj]] . 

The time evolution of the field operators is given by the form 

ip{x,t) = U+{t)il:{x,0)U{t) , 
expressed through the evolution operator U (t) satisfying the Schrodinger equation 

ij^U{t)^H[i;{x,0)]U{t). 

Respectively, the time dependence of the statistical operator p{t) stems from the Liouville 
equation, yielding 

m = u{t) m u+{t) . 

In any case, it is the grand Hamiltonian (13), which governs the temporal evolution of a 
nonequilibrium system. 



3 Broken Gauge Symmetry 

Now, let us specify the general notion of representative statistical ensembles, formulated 
above, for Bose systems, in which there exists the critical temperature Tc below which the 
global U{1) gauge symmetry becomes broken. For concreteness, we keep in mind a one- 
component system of spinless particles, characterized by the field operators satisfying the 
Bose commutation relations. 

Above the critical temperature Tc, the system is described by field operators ip = ip{r,t) 
and the conjugate ip"^ , being functions of spatial, r, and temporal, t, variables. These opera- 
tors generate the Fock space J-'{ip) on which all physical operators are defined. The related 
mathematical details can be found in the books [40-42]. The Hamiltonian energy operator 
H[il^] is a gauge-invariant functional of the field operators. The number-of-particle operator 
N[iIj] is normalized to the total number of particles =< N[ijj] >. Therefore the grand 
Hamiltonian is 

H[i;] = H[ij] - fiNiiJj] (T > Te) . 

This grand Hamiltonian, entering the statistical operator (6), characterizes the representative 
statistical ensemble for the normal Bose system, above the critical temperature, where the 
global gauge symmetry is preserved. 
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Below the critical temperature Tc, the global gauge symmetry becomes broken. The 
symmetry breaking is realized by the Bogolubov shift 

iP{r,t) iP{r,t) = r){r,t)+iP,{r,t) , (16) 

where r]{r,t) is the condensate wave function, while '0i(r, t) is the field operator of uncon- 
densed particles. The field operators V'l and ipl generate the Fock space which all 

physical operators are to be defined on. In Eq. (16), the condensate wave function ri{r,t), 
strictly speaking, is assumed to be factored with the unity operator Ijr in J-'{ipi). However 
here and in what follows, we shall use the common way of omitting the explicit appearance of 
the unity operator, in order not to make formulas too cumbersome. The condensate function 
ri{r,t) can also be termed as the coherent field, since the related coherent state I77 > is the 
vacuum state in the space This and other mathematical details can be found in Refs. 

[40,42,43]. 

Thus, instead of one field variable V'(r, t) above T^,, for the Rose system below Tc, where 
the gauge symmetry is broken, there arise two field variables, the condensate function (coher- 
ent field) 77(r, t) and the field operator ■0i(r, t) of uncondensed particles. These two dynamical 
variables are, of course, assumed to be linearly independent, being orthogonal to each other, 

lri*{r,t)M^,t)dr^O. (17) 

It is of principal importance to emphasize that the spaces for T > Tc and J-'{ipi) 

for T < Tc are mutually orthogonal [43,44]. The field operators ijj and ipi are defined on 
different spaces, J^{ip) and J-'{ipi), respectively, realizing two different unitary nonequivalent 
operator representations, with the Bose commutation relations [43,44]. As soon as the gauge 
symmetry is broken, one has to deal with the space J^{ijji). It would be principally incorrect 
to work, first, in the space J^{ip), accomplishing there some transformations, and then to 
pass to the space J^{ipi) by breaking the symmetry with the Bogolubov shift (16). As is 
shown in Ref. [6], such a procedure leads to wrong results. Prom the mathematical point 
of view, it is absolutely obvious that any manipulations must be accomplished in one given 
space, where all operators arc defined. 

Now, in the space Tiipi), we have two normalization conditions for two linearly inde- 
pendent field variables, rj{r,t) and V'i(r, t). The condensate function is normalized to the 
number of condensed particles 

No = J \v{r,t)\'dr, (18) 

which is assumed to be macroscopic, such that the condensate fraction Nq/N he nonzero in 
the thermodynamic limit. This normalization can be rewritten in the standard form of the 
statistical conditions (10) by using the operator 

^0 = Nol^ , (19) 

where Ijr is the unity operator in Then Eq. (18) transforms to the statistical 

condition 

No = < TVo > , (20) 

in which, as in what follows, the averaging is over the space J^(^i) . The second normalization 
is, clearly, for the number of uncondensed particles 

N, = < TVi > , (21) 
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with the corresponding operator 

Ni = J i/jl{r,t)Mr,t) dr . (22) 

The total number of particles 

N = <Ar> = No + Ni (23) 

is the average of the operator 

N = Jil^\r,t)i^{r,t) dr ^ No + Ni . (24) 

Normalization (23) follows from Eqs. (20) and (21). Therefore, among three normalization 
conditions, (20), (21), and (23), only two can be treated as independent. Generally, any 
combination of the pairs, for A^o and A^i, or for A^o and N, or for A^i and N, could be chosen. 
For the sake of symmetry, we prefer to choose the normalization conditions (20) and (21). 

When the gauge symmetry is broken, the average < ipi > may become nonzero. This, 
however, would mean that quantum numbers, as spin or momentum, are not conserved. In 
order to avoid this unpleasant situation, one has to impose an additional constraint 

< Vi(r,i) >= 0. (25) 

The latter can be reduced to the standard form of the statistical conditions (10) by defining 
a self-adjoint operator 

A[V^i] = J [X{r,t)4{r,t) + X*{r,t)M^,t)] dr , (26) 

in which A(r, t) is a complex function. Then constraint (25) can be rewritten as the quantum- 
number conservation condition 

< AfV'i] > = . (27) 

In this way, for the correct representation of a Bose system with broken gauge symmetry, 
we must work in the space and take into account three statistical conditions, (20), (21), 

and (27). The corresponding representative ensemble is constructed following the general 
procedure, formulated in Sec. 2. 

For an equilibrium system, according to Eq. (5), we have the information functional 

I[p] = IVp Inp + Ao (Trp - 1) + (3 (Trp H[ri, V'l] - e) 

- /5/xo (Trp No - No) - (Trp N^ - N^ - /?Trp A[V^i] . (28) 
Minimizing the latter yields the statistical operator 

p=i exp{-/3i7[r;,V^i]} , (29) 

with the partition function 

Z = IVexp{-/3//[77,Vi]} , 
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and the grand Hamiltonian 



H[r), ^i] = H[r), V^i] - iioNo - l^iN, - A[^i] , (30) 

in agreement with Eq. (7). 

For the general case of a nonequihbrium system, we have the Lagrangian 

L[77,Vi] = i[?7,V'i]-^[?7,V'i], 

in which 

Then, similarly to Eq. (12), the effective action is 

A[77,V'i] = / {Eh^l]-H][rJ,^^]} dt , (31) 

with the same grand Hamiltonian (30). Since we have now two linearly independent field 
variables, the extremization of action (31) gives two variational equations, for the condensate 
function, 

^%^ = 0, (32) 
Sr]*{r,t) ' ^ ' 

and for the field operators of uncondensed particles, 

5A[rj,ilJi] 



. (33) 



Equations (32) and (33), in view of the effective action (31), are equivalent to the evolution 
equations 

4r;(r,t) = ^3^ (34) 
dt ' ' Sr)*{r,t) ^ ' 

and, respectively. 

And, as usual, these equations are to be complimented by their Hermitian conjugate. 

Thus, the representative statistical ensemble {J-'{ipi), p,d/dt} for a Bose system with 
broken global gauge symmetry is defined in complete agreement with the general theory of 
Sec. 2. The dynamics and thermodynamics of such a system are governed by the grand 
Hamiltonian (30). It is only by accurately taking into account all imposed constraints that 
it is possible to correctly define the representative ensemble and to develop a self-consistent 
theory, avoiding any internal defects and paradoxes. The imposed statistical conditions (20), 
(21), and (27) lead to the grand Hamiltonian (30), with the Lagrange multipliers /xq, fJ^i, and 
X{r,t). The form of this Hamiltonian is more general than that of the trivial expression 
H — fiN. For a system with broken gauge symmetry, the number of the Lagrange multipliers 
in the form H — fiN is smaller than the number of imposed constraints. As a result, the 
problem becomes mathematically overdefined, which leads to the inconsistencies described 
in Sec. 1. 
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4 Self-Consistent Thermodynamic Relations 

For an equilibrium system, the condensate function does not depend on time, ri{r,t) = r]{r). 
The grand thermodynamic potential 

-T\nTrexp{-pH[7],ipi]) (36) 

is defined through the grand Hamiltonian (30). As is evident by this definition, 

dn on 

dfio djjii 
Varying potential (36) with respect to rjir) gives the equation 



5?7*(r) 5r]*{v) 

for the condensate function, in agreement with Eq. (34), when in equilibrium dr]{r)/dt — 
0. In order that the system with Bosc-Einstcin condensate be stable, the thermodynamic 
potential (36) is to be minimal with respect to the number of condensed particles iVo, so 
that 

5iVo - < ^ - ° ' ^^^^ 

which is the Bogolubov-Ginibre stability condition. Under normalization (18), the depen- 
dence on A^o comes through the condensate function r]{r), because of which 



01% 



Sn d7]{r) ^ Sn d7]*{r) 



5r]{r) ONo 5ri*{r) ONq 



dv . 



Therefore Eq. (38) is a direct consequence of Eq. (37). Of course, for a uniform system, 
Eqs. (37) and (38) identically coincide. 

The average densities of condensed, uncondensed, and all particles are 

_iVi _N 

respectively. The related fractions of condensed and uncondensed particles are 

From these definitions, one has 

P = Po + Pi , no + ni = 1 . (41) 

What one can usually fix in experiment is temperature T and the total density p. 
The system free energy writes as 

F = Q + poNo + piNi . (42) 

Recall that Nq = Nq(T,p) is such that to guarantee the system thermodynamic stability. 
From the equations of motion, together with the relation N — Nq + Ni, one finds A^i = 
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Ni{T,p). As is mentioned above, in experiment, only the total number of particles can be 
fixed, except temperature and volume. Hence, the free energy can be written as 

F^n + liN , (43) 

which can be considered as the definition of the system chemical potential /j,. Comparing 
Eqs. (42) and (43) gives 

II = iiono + jJ-iUi . (44) 

Then one can define the free energy F — F{T,V,N) as a function of temperature, volume, 
and the total number of particles, with the differential 

dF ^ -S dT - P dV + II dN , (45) 

in which 5* is entropy and P is pressure. From here, all measurable thermodynamic quantities 
are calculated in the standard way. For example, the system chemical potential is 

As is evident, Eqs. (44) and (46) are identical l)y the above definitions. 

At zero temperature, the free energy F — E — TS reduces to the internal energy 

E=< H[r], V'l] > = < H[ri, V'l] > +iiN . (47) 

Therefore the chemical potential (46) becomes 

The grand potential (36) is a function Q = Q(T', V, jj) of temperature, volume, and 
chemical potential, with the differential 

dVL^ -S dT - P dV - N dix . (48) 

It is easy to check that the thermodynamic and Gibbs entropies coincide, so that 

5 = -|^ = -IVplnp. (49) 

This immediately follows from the direct differentiation of the grand potential (36) and the 
form of the statistical operator (29). 

It is also possible to prove (see the proof in the Appendix A) , that the standard relation 

f = /?A^W (50) 

holds, in which A^(iV) is the dispersion of the number-of-particle operator N — Nq -\- N^. 
The dispersion of a self-adjoint operator A is 

A2(i) = < > - < i >2 . 

In this way, all thermodynamics quantities and relations are defined self-consistently. 
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5 Operator Equations of Motion 



The equations of motion for the field variables ?7(r, t) and ipi{r, t) are given by Eqs. (34) and 
(35). To specify them, let us take the Hamiltonian energy operator in the standard form 



^[77, Vi] = / ^^^(r) (- ^ + ^) dr 

+ 1 J V'Hr)^^(r')*(r - r')^(r')^(r) drdr' , (51) 

in which '?/'(r) = V'(r, t) is the shifted field operator defined in Eq. (16). In order to avoid 
cumbersome notations, we omit the explicit dependence on time of the variables rj{r) = rj{r, t) 
and ^i(r) = ipi{r, t), though this dependence is assumed. In what follows, such an abreviated 
notation will be often used, where it does not lead to confusion. The external potential 
U — U{r,t) can, in general, depend on time. The interaction potential $(r) is arbitrary, 
provided it is symmetric, such that $(r) = $(— r), and enjoys the Fourier transformation. 
The grand Hamiltonian (30) can be written as a sum of five terms. 



(52) 



n=0 



whose order is labelled by the number of operators ipi in their products. The zero-order 
term contains no ipi, being 



^ Jn*{r)(^-^ + U- /.o) dr 
+ y - r'Mv')\'\n{r)\' drdr' . 



(53) 



To satisfy the quantum-number conservation condition (27), the Hamiltonian must not con- 
tain the terms linear in ipi. This prescribes to take the Lagrange multipliers in Eq. (26) 
as 

^\u + j $(r-r')|77(r',t)rrfr' 



\{r,t) 



2m 



(54) 



so that H^^'' = 0. The necessity of removing the terms linear in ijji or V^j, by means of condi- 
tion (54), in order to satisfy constraint (25), can be easily proved for quadratic Hamiltonians, 
involving linear terms, by employing the method of canonical transformations [42,45]. In 
the general case of arbitrary Hamiltonians, the proof of Eq. (54), yielding the cancellation 
of linear terms, is presented in the Appendix B. 
The second-order term is 



2m 



+ U - Hi] 'ipi (r) dr 



+ 



$(r - r') f|r](r)|2^I(r')t/'i(r') + r/*(r)r;(r')^/'I(r')^i(r) 

drdr' . 



+ \ rr{r)rf {r')Ur')Ur) + \ r){r)r){r')4{r')4{r) 



(55) 
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Respectively, we have the third-order term 

^(3)^ J $(r-r') [r;*(r)V^l(r')V'i(r')V^i(r) +V^l(r)V^l(r')V^i(rXr)] drdr' 
and the fourth-order term 

= ^ y V'I(r)V^l(r')$(r - r>i(r')V'i(r) drdr' . 
Prom Eq. (34), we obtain 

t ^ ry(r, t)^(^-^ + U-i,o^ 77(r, t) + J $(r - r') [\v{r')\'ri{r) + ^(r, r') 

where the last term in the integrand is the correlation operator 

X{r,r')=4{r')M^')rj{r) 

+ 4{r')7j{r')ij,{r) + rj* {r')Mr')M^) + iplir')Mr')Mr) ■ 
And Eq. (35) yields 

d / \ 

^ — ^i(r,i)= i-— + U-i^,Ui{r,t) 



dt 



I 



+ / $(r - r' 



r^(r')|Vi(r)+r^*(r')ry(r)V'i(r') + r^(r')r^(r)Vl(r')+X(r,r') dr' . (60) 



(56) 



(57) 



dr' , (58) 



(59) 



Again, for brevity, the time dependence is not explicitly shown in the integrals of Eqs. (59) 
and (60). Also, recall that in the operator equation (58), the condensate function is assumed 
to be factored with Ijf, the unity operator in 



6 Local Conservation Laws 

The equations of motion (58) and (60) are derived from the variational principle of action 
extremization. Therefore they guarantee the validity of all local conservation laws on the 
operator level. As an illustration, let us consider the time variation of the local densities. 
The local condensate density is 

Po{r,t) = \rj{r,t)\' , (61) 

and the local condensate current is 

jo(r, t) = -^ [r?*(r, t) Vrj{r, t) - rj{r, t) Vr^*(r, t)] . (62) 

Respectively, we define the operator density of uncondensed particles 

p,{r,t)=4{r,t)M^,t) (63) 

and the related current operator 

i,{r,t) = -^{4{r,t) VM^,t) - [v4{r,t)]M^,t)} . (64) 
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Differentiating the condensate density (61), we find 

^ Po(r, t) + V- jo(r, t) = -f(r, t) , (65) 

where the source operator is 

f(r, t) = J $(r - r') [^(r, r') + ^+(r, r')] dr' , (66) 

in which 

i?(r, r') = ir]*{r) [iPl{r')r]{r') + 77*(r>i(r') + ^I(r')(/'i(r')] Vi(r) ■ 
And the time variation of density (63) yields 

^Pi(r,t)+V-ji(r,t)=f(r,i). (67) 
For the total operator density 

p{r,t) = po{r,t) + pi{r,t) (68) 

and the total operator of current 

j(r,t)=jo(r,t)+ji(r,t), (69) 
we obtain the continuity equation 

^p(r,t) + V-j(r,t)=0. (70) 

In the same way, one can derive any other local conservation laws, following the standard 
procedure [40,45] and employing the equations of motion (58) and (60). 

7 Condensate Wave Function 

The equation for the condensate wave function r]{r,t) follows from averaging Eq. (58) over 
J-'{ipi). To this end, we need to introduce several notations. The normal density matrix is 

pi(r,r') = < Vl(r')V'i(r) > ■ (71) 

Under the broken gauge symmetry, there appears the so-called anomalous density matrix 

(7i(r,r') = < V'i(r')V'i(r) > ■ (72) 

The density of condensed particles is po(r), defined in Eq. (61), and the density of uncon- 
densed particles is 

pi(r) = pi(r,r) = < Vl(r)V^i(r) > . (73) 
We shall also need the anomalous average 

(71 (r) = ai(r,r) = < ^i(r)^i(r) > . (74) 
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The absolute value |(Ti(r)| has the meaning of the density of pair-correlated particles [43]. 
The total density 

p(r) =po(r)+Pi(r) (75) 

is the average of the operator density (68). According to the normalization conditions (20) 
and (21), the partial densities are normalized to the number of corresponding particles, 

No^l po(r) dr, N^^ j pi(r) dv . (76) 
One more notation is the triple correlator 

e(r,rO = < Vl(r')Vi(r')V'i(r) > ■ (77) 
Using these definitions, the average of the correlation operator (59) becomes 

< X(r,r') > = pi(r')77(r) +pi(r,r')77(r') + ai(r,r')77*(r') +e(r,r') . 
Finally, averaging Eq. (58) yields the evolution equation for the condensate wave function 

d / \ 

+ j $(r-r')[p(r')r/(r) + pi(r,r')r/(r') + ai(r,r')r/*(r')+e(r,r')] dv' . (78) 

This is the general equation for an arbitrary Bose-condensed system. This equation can also 
be represented in a shorter, though not explicit, form 

I ^ r7(r, t)^[-^ + U-ii^ rj{r, t) + J $(r - r') < ^i;^(r')i^{r')^(r) > dr' , 

in which one has to substitute the shifted operator ip — rj + ipi. 

We may notice that Eq. (78) is invariant under the transformation 

r/(r,t) ^r/(r,t)e*"*, ^i(r,t) ^ ^,{r,t)e''^\ /io ^ /io + a • (79) 

This means that there exists a freedom in choosing the phase factor exp(iQ;t) of the conden- 
sate function. But the phase factor becomes fixed by defining the Lagrange multiplier /xq. 
The arbitrariness in the condensate phase imphes that the stationary solutions for r]{r,t) 
would be proportional to an undefined factor e'"*. By fixing the multiplier fj,Q, we require 
that in equilibrium the condensate function would not be dependent on time, that is, 

d 

— r]{r) = (equilibrium) . (80) 
In equilibrium, Eq. (78), according to condition (80), reduces to the eigenvalue problem 



/^o?7(r) 



r){r) 



+ 



J $(r - r') [pir'Uv) + p^{r, v')r]{v') + a^{v, r')r]*{r') + ^(r, r')] dv' , (81) 
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defining, together with the normahzation condition (76), the eigenfunction ?7(r) and the 
eigenvalue ^q. Let us emphasize that without the normahzation condition (76) the con- 
densate function cannot be uniquely defined from Eq. (81). For example, if the system is 
uniform, with no external potential U ^ 0, when ?7(r) = r], p(r) = p, and r] = then Eq. 
(81) gives 

e(r,0)- 



/io = p$o + j $(r) 

where $o = / ^{r)dr, or in the compact form 

1 



pi(r,0) + (7i(r,0) + 



dr , 



(82) 



1^0 



$(r) < VUO)V'(0)V'(r) > dr . 



Expression (82) is valid for an arbitrary uniform equilibrium system with an interaction 
potential $(r). No approximations have been used in obtaining Eq. (82). 

In recent days, the physics of dilute Bose gases is intensively explored (see the book [46] 
and review articles [31,47-49]). The interaction potential for these gases is taken in the 
contact form 



$0 = 47r — , 
m 



(83) 



where is the scattering length. With this interaction potential, the eigenvalue problem 
(81) for a nonuniform system becomes 



/^o^(r) 



2m 



+ U{v) 



r){r) 



+ $0 {[p(r) + Pi(r)] r] + a,{r)rj*{r) + ^(r, r)} . (84) 

For a uniform system, /iq is given by Eq. (82), which, in the case of the contact potential 
(83), yields 



p + Pl + (7i + 



(85) 



where (7i = cri(r, r) and ^ = ^(r, r). 

When one assumes that the Bose gas is dilute, being characterized by the contact po- 
tential (83), the temperature is zero, and the interaction is so weak that the condensate 
depletion can be neglected, so that all particles are condensed, then Nq = N, Ni = 0, and 
Pi = 0"! = ^ = 0. In this approximation, because of relation (44), the multiplier fiQ = fi 
coincides with the system chemical potential. The eigenvalue problem (84) reduces then to 
the Gross-Pitaevskii equation 



— + C/(r) + $okr)r 



77(r) = iir){r) 



But when the condensate depletion is not negligible, then //q, according to Eq. (44), is not 
the same as /i. 
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8 Uniform Bose system 

Let us illustrate in more detail the application of the representative ensemble with the grand 
Hamiltonian (52), to an equilibrium uniform system. Then the field operators of uncondensed 
particles can be expanded in plane waves, 

^i(r) = l^afc<^fc(r) , (86) 

where v^A;(r) = e^^'''/\/V. The condensate function reduces to the constant r/ = ^/po, while 
the condensate multiplier /iq is given by Eq. (82). 

In the momentum representation, the momentum distribution of particles 

nk = < alttk > (87) 

is usually termed the normal average, as compared to the anomalous average 

^k = < a-ka-k > ■ (88) 

The normal and anomalous density matrices (71) and (72) take the form of the expansions 

r') = nkV^k{r)vl{r') , cri{r, r') = ^ akipk{r)cpl{r') , (89) 

k^Q fc^O 

in which the properties 

are taken into account. The diagonal elements of the matrices in Eq. (89) give the densities 
Pi = Pi(r, r)^^ Y^Uk, (Ti = (Ti(r, r) = ^ (jfe . (90) 

^ k^O ^ fe^O 

The interaction potential is assumed to allow the Fourier transformation 

^(r) = 7 E ^fee^''"' ^fe = / ^(Oe-'*^"- dr . (91) 

k 

Then we find the following terms of the grand Hamiltonian (52). The zero-order term (53) 
becomes 

= Qpo^o-/^o)a^o. (92) 
The first-order term H^^ = is automatically zero. The second-order term (55) is 



r ^2 "I 
^^^^ = 11 ^ + Po(^o + ^fe) - A*i 4afc + - ^ po^fe (4«-fc + ■ (93) 



For the third-order term (56), we have 



2 

^ k^o 



H^^^ ^ ]lyY! % (alak+pa-p + alpul^pak) , (94) 

kp 
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where the prime on the summation symbol imphes that k 7^ 0, p 0, k + p 7^ 0. The 
fourth-order term (57) is 

^^^^ = ^577 EE' $5«I«J«P+9«fc-<z ' (95) 

q kp 

where the prime shows that k^^O, p^^O, p + qT^O, k— q^^O. 

To be more specific, it is necessary to resort to some approximation. The natural mean- 
field approximation for a system with broken gauge symmetry is the HFB approximation. 
The latter is usually blamed to display an unphysical gap in the spectrum, because of which 
it is qualified as gapful (see detailed discussion in Refs. [7,31,32]). However, as is explained in 
the Introduction, this defect comes into play only owing to the usage of a nonrepresentative 
statistical ensemble. But for the representative ensemble, with the grand Hamiltonian (52), 
there appear no such problems. Below, we show this for the case of arbitrary temperature 
and interaction potential $(r). 

We apply, in the standard way [40,45], the HFB approximation to the third- and fourth- 
order products of the operators a^. Then the third-order term (94) becomes identically zero, 
because of the condition < >= 0, 

H^^'^ = . (96) 
And for the fourth-order term (95), we get 

k^o ^ ^ 



1 r 1 1 

+ E ^k+pdldp + 2 {cTk+pala-p + (^k+pO'-pO'p) - 2 i^k+pUp + (Jk+p(Tp) 

k,p^O 

It is convenient to introduce the notations 



(97) 



u;k = — + p^o + Po^k + TT E ^p^k+p - Pi (98) 

^ p^O 

and ^ 

Afe = po^k + ^7 E ■ (99) 

Since the interaction potential $(r) = $(— r) is symmetric and real, we have the properties 
(7fe = (7^ = (7_fe and Afc = A* = A.^.. 

Summing up all terms in the grand Hamiltonian (52), we obtain in the HFB approxima- 
tion ^ 

Hhfb = Ehfb + E + 9 E + «-fe«fc) > (100) 



2 

k^O ^ kjtO 



where the nonoperator term is 



Ehfb = H^''^ - ^pI%V - ^ E ^k+pirikUp + akap) . (101) 
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The quadratic form (100) can be diagonalized by the Bogolubov canonical transformation 
[3,4,45] 

flfc = Ukbk + v*_i^bl^ . 

As a result, the grand Hamiltonian (100) is transformed to the Bogolubov representation 

HB^EB + Y^Skblbk, (102) 

in which 

Eb = Ehfb + I - ^k) . (103) 

For the operators bk, one has the properties 

< 6fe > = < > = , <blbp>^ SkpTTk , (104) 
with the momentum distribution of quasiparticles 

TT, = < blbk > = (e^^'= - l)"' = ^ Mh ■ 



(105) 



The coefficient functions Uk and v^, and the spectrum are defined by the Bogolubov - de 
Gennes equations 

{(jjk - Sk)uk + AkVk = , AkUk + (wfc + £fc)vfe = , 
with the normalization condition — v^. — 1. This gives 

^2 _ '^k + ^2 _ ^ ~ 

And the Bogolubov spectrum is 



Sk = ^col - Al . (106) 



The existence of the Bose-Einstein condensate, as is known, requires that the spectrum 
(106) be gapless, such that 

lim£jfc = 0, (107) 

k—*0 

with the stability condition Sk > 0. Then Eq. (106) yields 

III ^p% + ^ ^(np - ap)% . (108) 

As will be shown in Sec. 10, this value is in exact agreement with the Hugenholtz-Pines 
relation. This should be compared with the condensate potential (82), which in the HFB 
approximation, when ^(r,r') = 0, becomes 

Ho ^p^o + ^ Y.{np + ap)% . (109) 

In the particular case of the contact potential (83), we get 

Ail = (P + Pi - f^i)*o (110) 
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and, respectively, 

/xo = (p + Pi + (7i)$o ■ (111) 

In any case, generally, //q 7^ l^i- The multipliers /Xq and Hi become equal only in the 
limit of the asymptotically small condensate depletion, which corresponds to the Bogolubov 
approximation [1,2], when both pi and ui in Eqs. (108) to (111) are neglected. Then, 
clearly, Ho = Hi = po^o, and in view of relation (44), /iq = pi = fi. However, as soon as the 
condensate depletion is not negligible, //o 7^ A*i- 

With the multipher (108), Eq. (98) can be represented as 

c^fc = ^ + Po^k + 77 J2 {np^k+p - np% + (7p$p) . (112) 

For the long- wave limit of A^, we introduce the notation 

A = limAfc = mc^. (113) 

fc— >o 

Prom Eq. (99) it follows 

A = po^o + 77 ^p^P ■ (114) 

Then the long-wave limit of spectrum (106) is explicitly of the phonon type, e ~ ck, as 
A; — > 0. According to the notation (113), the sound velocity is 

(115) 

Strictly speaking, the sound velocity is defined as c = ^jA/m*, with an effective mass to be 
given in Sec. 10. For short-range interactions, m* = m. Using Eqs. (113) and (114), we 
may rewrite Eq. (112) in the form 

uJk^ — +mc^ + poi^k - *o) + J np{^k+p - %) ■ (116) 

For the particle momentum distribution (87), we find 

n, = ^coth(^) - -, (117) 
2£fe \2TJ 2 ' ^ ' 

while for the anomalous average (88), we obtain 

Analyzing the behaviour of and ak as functions of momentum k, we find [15,32] that 
|(Tfc| ~ nfc for A; — > 0, while \ak\ ^ Uk for large k. Therefore in no sense the anomalous 
average ak can be neglected, as compared to the normal average n^. 

The derived equations simplify for the contact potential (83). Then = ^k — 

A = (po + (Ti)$o , 
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and spectrum (106) takes the classical Bogolubov form 



A;2 Y 



^W+^-j . (119) 

The grand potential (36) in the HFB approximation is 

n^Es + TV I ln{l-e-^^^) (120) 

From this and other equations, obtained above, all thermodynamic characteristics can be 
calculated. 

9 Condensate and Superfluid Fractions 

We shall concentrate our attention on the most interesting characteristics of the Bose- 
condensed system, on the condensate and superfluid fractions. The condensate fraction 

no = 1 - ^ (121) 
P 



can be found by calculating the density pi of uncondensed particles, given by the integral 

dk 
(2^ 



P^-jn.j^,. (122) 



To define the superfluid fraction, one considers the reaction of the system to the boost 
with a velocity v. In the laboratory frame, the field operators of a moving system are 
represented by means of the Galilean transformation 



■^'^(r, t) = ■?/'(r — vt, t) exp |i ^mv • r — J ■ (123) 

The Hamiltonian — H^^], in terms of the new field operators (123), is expressed through 
the Hamiltonian H — H[^l}\, in terms of the old operators ■?/;, as 

H,^H + I V'^(r) ^-iv ■ V + ^ j ^(r) dr . (124) 

The total system momentum operator becomes 

= P + Nmv , (125) 

where 

P = I i,^{r){-iV)ip{r) dr (126) 
is the old momentum operator. The average total momentum of the moving system is 

< P^ >^ = Trp„P, , (127) 
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with the statistical operator 



^ ^I-^ • (128) 



The superfluid fraction is defined as 

1 d 

It is possible to show (sec, e.g., the detailed derivation in Refs. [47,49]) that the superfluid 
fraction (129) for an arbitrary system can be represented in the form 

Us^l (130) 



in which 



Q = =-^ (131) 
^ 2mN ^ ' 



is the dispersed heat, with the momentum dispersion 

A2(P) = < p2 > - < P 

where the averages are with respect to the grand Hamiltonian ['?/']. 

Formula (130) for the superfluid fraction is valid for any system, equilibrium or nonequi- 
librium, uniform or not. For a system in equilibrium, the total momentum is zero, < P >= 0, 
so that the dispersed heat (131) is 

< P2 > 

Q = . ■ (132) 

For a uniform system, 

< P' > = E (k ■ P) < > , (133) 

kp 

where fik = ajtOfe- In the HFB approximation, 

< fikfip > = rikUp + SkpUkil + Uk) + S_kpcrl . (134) 
Then the dispersed heat (132) becomes 

Substituting here expressions (117) and (118), and using the equality coshx—1 = 2sinh^(x/2), 
so that 

1 



Uk + nl- al 



4sinh2(/3£fe/2) ' 
we get from Eq. (135) 

p_ 1 r ^'^^ n^fi^ 

(47r)2mp Jo smh'{(3sk/2) ' ^ ^ 
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Note again the importance of the anomalous average au- If one would omit the latter in the 
dispersed heat (135), one would get a senseless divergent quantity, while accurately taking 
account of results in the well-defined convergent integral (136). 

To demonstrate explicitly the behaviour of the condensate and superfluid fractions, let 
us resort to the contact potential (83). Then Eq. (99) yields 



and Eq. (116) gives 



mc 



2m 



Equation (115), defining the sound velocity, can be written as 

mc^ = (po + (7i)$o , 



where 



0-1 



dk 



The density of uncondensed particles (122) can be represented in the form 



Pi = 



(mc)' 
37r2 



1 + 



3 



2a/2 Jo 



(VTT 



1/2 



, / mc ' 
coth — — X 
\ 2T ; 



dx 



(137) 



(138) 



(139) 



(140) 



(141) 



which is a well-defined convergent integral. 

The anomalous average (140) can be written as a sum of two parts 



(Ti — ao — 



mc 



2s, 



coth(^)-l 



dk 



in which 



^ _^ /• 1 dk 

J 2^kWy 



(142) 



(143) 



The second integral in Eq. (142) is convergent. But the integral in Eq. (143) ultravioletly 
diverges. However, this divergence is known to be unphysical, being simply caused by the 
usage of the contact interaction potential. This and other similar divergences could be 
removed by employing more realistic interaction potentials. For example, a common choice 
is a Gaussian type potential [50]. Another known way of removing such divergences is 
by using the analytic regularization in one of its variants, such as the subtraction scheme, 
zeta regularization, or dimensional regularization (see details in Ref. [31]). This kind of 
regularization is asymptotically exact for weak interactions and is universal in the sense that 
it applies to any short-range potential with a scattering length [31]. Using the dimensional 
regularization in its region of validity, we obtain 



1 dk 

2i^ 



^ Jmpo^o 



Then Eq. (143) results in 



[mc) 



^Jmpo^o . 



(144) 
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The anomalous average (142) can be represented as 



(7i = (7o 



2^2 7r2 Jo y/T 



coth 



mc 



X 



dx 



(145) 



And for the dispersed heat (136), we obtain 

r~ {^/TT^ - 1)^/^ xdx 



Q 



[mc] 



1 

Jo 



A/2(27r)2mp Jo y/l + x^ smh.'^{mc'^x/2T) 
At low temperatures, such that 

T 

«1, 



(146) 



equations (141), (145), and (146) yield 

^3 {mcf f T 



mc 



12 V 



Q 



(mc)^ ( T 



(7i ~ (7o 



12 \mc 



7r2(mc)5 / T 



15mp \m(? 
Therefore, the condensate fraction (121) behaves as 



no ~ 1 



{mcf {mcf ( T 



Sn^p 12p \mc^ 
And the superfluid fraction (130) becomes 

27r\mcf f T 
n., ~ 1 - ^ ' I 



(147) 



45p Kmc 



(148) 



When the interaction is weak, or the condensate fraction uq tends to zero at the critical 
temperature Tc, as a results of which, c — > 0, so that 



mc 



«1, 



then Eqs. (141), (145), and (146) lead to 



rp . 3/2 



Pi ~ P — + 



(mcj"" 
"3^ 



m^cT 



CTi ~ CTq 



2tt 



with the critical temperature 



mc 

7;J ~ c(3/2) Vt;; 'tT 



2tt 



m 



P 

.am. 



2/3 
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coinciding with that for the ideal Bose gas, as it should be in the case of the mean-field 
approximation with the contact potential (83). Here ('(•) is the Riemann zeta function. 
Then the condensate fraction tends to zero at Tc as 



together with the superfluid fraction 

Both, no and n^, tend to zero from the left, demonstrating the second-order phase transition 
at T(.. The point why in the present case T^. has to be the same as the ideal Bose gas critical 
temperature is explained in the Appendix C. 



10 Green Function Equations 

In order to conclude the proof of the complete self-consistency of the developed approach, 
based on the introduced representative ensemble, we need to show that for an arbitrary Bose 
system there exist the Green-function equations of the usual form and that the Hugenholtz- 
Pines relation follows for any uniform system. 

The first-order Green function is defined in the standard way [3,4,7,9] as the matrix 
^(12) = [Ga/3(12)] with the elements 

G'n(12) = < ^^^1(1)^1(2) > , 6-12(12) ^ -i< f ^i(l)V^i(2) > , 

6-21(12) = -i< f^I(l)^I(2) > , 622(12) = -i< f^I(l)^i(2) > , (151) 

in which the common abbreviation is employed, denoting the set {rj, t^} by a single number 
J, and T is the chronological operator. 

For what follows, we shall need the notation for the triple field operator 

^(123) = V^i(l)V^i(2)^I(3) + r;(l)V^i(2)^I(3) + V^i(l)r;(2)V^I(3) 

+ V'i(l)V'i(2)77*(3) + 77(1)77(2)V'I(3) + r^{l)^r{2)r^* {?,) + Vi(l)77(2)77*(3) ■ (152) 
One may notice that 

^(123) = V^(l)V^(2)V^t(3) _ ^(i)^(2)r;*(3) , 

where = 77 -|- -01 is the shifted field operator. However, for practical application, we need 
the explicit form (152). 

The second-order Green function is a matrix 5(1234) — [i?Q,/3(1234)] with the elements 

Sn(1234) = - < f*(123)^/'l(4) > , Si2(1234) = - < f *(123)V'i(4) > , 
^21(1234) = - < r*t(123)^I(4) > , ^22(1234) = - < r*^(123)V^i(4) > . (153) 
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For the general form of a retarded interaction potential 

$(12) = $(ri-r2)5(ti2 + 0), 
in which tu — h — t2, the self-energy is introduced through the equation 

$(13)5(1332) d{3) . 



S(13)G(32) d{3) = 
The equation for G{12), by defining the inverse propagator 

G-^(12) U{1) + pi,^ 5(12) 1 - E(12) , 



(154) 



(155) 



(156) 



in which 



" 1 


" 




" 1 


" 




i = 







-1 







1 



can be represented in the form 



J G-\13)G{32) d{3)^S{12) 



1 . 



(157) 



To solve Eq. (157), one may invoke perturbation theory, starting with an available 
approximate Green function Gapp, corresponding to an approximate E^pp, such that 

J G-^^{13)Gappi32) d{3) = 6(12)1 . (158) 

Then Eqs. (157) and (158) can be transformed to the Dyson representation 

Gil2) = Gapp{12) + J Gappil3) [E(34) - E„pp(34)] (7(42) d{3A) , (159) 

which is convenient for using perturbation theory. 

The above equations for the Green functions are valid for an arbitrary nonequilibrium 
and nonuniform Bose system. If the considered system is equilibrium and uniform, one 
passes to the Fourier transforms ^(k, a;) and E(k, a;), employing the symmetry properties 



Gi,{k,-uj)=G22{k,uj) , Gi2(k,-a;) =G'2i(k,cu) =G'i2(k,w) , 
G'a/3(-k,u;) = GaoiK^) 

and, respectively, 

Eii(k, -uj) = E22(k,a;) , Ei2(k, -a;) = E2i(k,a;) = Ei2(k,a;) 

Ec,/3(-k,u;) = Ea/3(k,u;) . 

A detailed discussion of these relations can be found in Ref . [4] . 
Let us introduce the notation 



(160) 



(161) 



D{k,uj) = El,{k,uj)-G^lik,u;)G^,\k, -co) , 



(162) 
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in which 



Gil{k,u;)^u; - T,n{k,u;) + Hi . 

2m 



Then the solution to Eq. (157) can be written as 

a; + /i;V2m + Eii(k,u;) -//I E^iKu) 
G'ii(k,a;) = — ^ , Gi2{k,uj) ^ - —— -, (163) 



D{k,u;) 



D{k,uj) ' 



with the denominator (162). 

For the Green functions Ga/sik, 0) at zero energy there is the Bogolubov theorem [4] 
rigorously proving the validity of the inequalities 



|Gii(k,0)| > 



mriQ 



|Gn(k,0)-Gi2(k,0)| > 



From Eq. (163), we have 



Gii(k,0) 



where 



_ A;V2m + Sii(k,0) - /^i 
" i?(k,0) 

L'(k,0) = E?2(k,0) 



^12(1^,0) 



Si2(k,0) 
L'(k,0) 



n 2 



2m 



+ Eii(k,0) -/ii 



Consequently, 

6-11(^0) -6-12(^0) = 
Using this in Eq. (165), we get 



2m 



+ Ei2(k,0)-En(k, 0) 



-1 



1^1 - — + Ei2(k,0)-Sn(k,0) 
zm 



< 



mno 



Setting in the latter equation k — > 0, we come to the Hugenholtz-Pines relation 

//i = En(0,0)-Ei2(0,0) . 
Note that in the HFB approximation we have 

En(0, 0) = (p + po)^o + 7 E ^p^p ' ^12(0, 0) = po^o + ^ E ■ 



V 



(164) 
(165) 



(166) 



Therefore Eq. (166) gives exactly the form of //i in Eq. (108). 

Relation (166) guarantees that the system spectrum is gaplcss and of the phonon char- 
acter. The spectrum is given by the zeros of the Green functions (163), that is, by the 
equation 

D{k,ek) = 0, (167) 
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with Eq. (162). Equation (167) can be rewritten as 



£fe = ^ [En(k, £fe) - E22(k, s^)] + .Jjf^Y^Jy^) , (168) 

where the notation 

t^fc = ^ + 2 Pii(k, ^k) + S22(k, £fe)] - /ii 

is used. From Eq. (168), it follows that Sk ^ 0, as k ^ 0. Moreover, when the sys- 
tem is isotropic, its self-energy SQ^(k, e^) depends only on the scalar k"^. This implies the 
asymptotic, as /c — > 0, expansion 

E,;3(k,£ife)~E,^(0,0) + E^^A;% 

in which 

d 

^'"/5^l™0 dk^ ^«/3(k,£fc) . 

Using this expansion in spectrum (168) gives 



Sk^ck, c=W^Ei2(0,0), (169) 
V m* 

which is the phonon spectrum, with the sound velocity c, where the effective mass is 

" 1 + m(E'n + S'22 - 2E'i2) ' 

It is useful to compare the Lagrange multipliers and given by their general ex- 
pressions (82) and (166), which are valid for an arbitrary equilibrium uniform Bosc system. 
These expressions are exact, with no approximations being involved. As is seen, anomalous 
averages enter hq with the sign plus, while the anomalous self-energy enters /ii with the 
sign minus. This is why, in general, no cannot coincide with /xi. The HFB forms (108) 
and (109), or (110) and (111), are particular illustrations. It is easy to check that yUg co- 
incides with Hi solely in the Bogolubov limit of asymptotically weak interactions, when 
Eii(0, 0) — >• (p-|- po)^0) ^12(0, 0) — s> po^O; and iiq ^ Hi ^ 11 = p$o- But in any higher-order 
approximation, //q 7^ A*i- The assumption that /Iq would be the same as /ii would make the 
theory not self-consistent and would return us back to the Hohenberg-Martin dilemma of 
conserving versus gapless theories. 



11 Discussion 

The main message of this paper is the necessity of employing representative ensembles for 
correctly describing statistical systems. A representative ensemble takes into account all 
imposed constraints and additional conditions that allow for a unique description of the con- 
sidered system. It is only using a representative ensemble makes the theory self-consistent. 

In the Bose system with broken global gauge symmetry, realized by the Bogolubov shift, 
there are two particle components, corresponding to condensed and uncondensed particles, 
with two related normalization conditions for A^o and A^i. This requires to introduce two 
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Lagrange multipliers, jio and //i, which makes the theory completely self-consistent in any 
approximation. 

It is worth recalling that the introduction of several Lagrange multipliers is rather com- 
mon for spin systems. There, the order parameter is the average spin, which, generally, is 
a three-component vector. The role of the effective chemical potential for spin systems is 
played, as is well known, by an external magnetic field, which is also a three-component 
vector. Hence, the number of effective chemical potentials for spin systems is equal to the 
number of components in the order parameter. Only then one is able to unambiguously 
define the average spin. 

The suggested approach, introducing two Lagrange multipliers, does not contradict our 
physical understanding that the standard experiments fix, as independent variables, tem- 
perature, volume, and the total number of particles. To emphasize this once again, let 
us turn to the definition of the grand thermodynamic potential (36), from which it fol- 
lows that it is a function fl = fl{T,V, /iq, /ii), so that the free energy (42) is a function 
F — F(T,V, No, Ni). The Lagrange multipher //q is defined from the stability condi- 
tion, yielding //q = fXo{T,V, Nq, Ni), according to Eqs. (81), (82), (84), (85), (109), and 
(111). The Lagrange multipler satisfies the Hugenholtz-Pines theorem, which gives 
f^i = fJ'iiTjV, No, Ni), in agreement with Eqs. (108), (110), and (166). The number of 
uncondensed particles is found from the direct calculation of the average (21) expressed 
through Eqs. (90), (122), (141), and like that, resulting in N^ = Ni{T,V, No, N). Prom 
here, since A^o = ^ - ^^i, it follows that A^o = No(T, V, N). Substituting it back to A^i, one 
has A^i = Ni(T, V, N). Using these in the expressions for and fii, one gets = /xo(T, V, N) 
and Hi = fj,i(T,V, N). Then, from Eq. (44), it is evident that fj, = fi(T,V,N), or, inverting 
the latter relation, one has N = N{T, V,/x). Using this, the Lagrange multipliers /^o and /Xi 
can be expressed as functions fj,o — l^o{T, V, /i) and /xi = /ii(T, V,fi). Substituting this into 
the grand potential, we have Q — Q{T,V,ijl), in line with Eq. (48). Respectively, the free 
energy becomes a function F = F{T, V, N), in accordance with Eq. (45). 

Thus, at the end, we work with the standard variables T, V, and A^, which are usually 
fixed in experiments. All observable quantities are also expressed through these variables. 
So that the suggested approach is absolutely self-consistent, mathematically correct, and in 
agreement with physics. 

Several words are to be said with regard to the phase-transition order of Bose-Einstein 
condensation. This transition is known to be of second- order, which is firmly based on sev- 
eral facts. First, there exists a general explanation, independent of the coupling strength, 
that this transition is of second-order. This can be found in the book by Patashinsky and 
Pokrovsky [51] (Chapter X, Section 2). As is also well known, the Hamiltonian of Bose sys- 
tems is mathematically equivalent to what in quantum theory is termed the c/?^ model. The 
phase transition in this model has been studied in numerous works using the renormaliza- 
tion group approach, exhibiting the second-order transition [52]. The superfiuid transition 
in liquid ^He, which is believed to be accompanied by the Bose-Einstein condensation, is 
also a second-order transition. A large body of experimental data on measuring the con- 
tinuous temperature dependence of the condensate fraction in superfiuid helium has been 
summarized by Wirth and Hallock [53]. There exists abundant literature, both theoret- 
ical and experimental, on Bose-Einstein condensation in dilute trapped gases (see review 
works [46-49]) and references therein) and there are several computer simulations of this 
process [54-56] . Though Bose condensation in traps is smeared out by finite-size effects, the 
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subsequent increase of the number of particles unambiguously demonstrates that the con- 
densation approaches the standard second-order phase transition. In addition, there have 
been many Monte Carlo calculations for uniform Bosc systems with various interaction po- 
tentials. All these calculations, summarized in the review articles [57-59] , clearly prove that 
Bose-Einstein condensation is a second-order transition. So, the second order of the Bose 
condensate transition has been established without any doubt. This especially concerns the 
general theoretical investigations [51,52] and rigorous Monte Carlo calculations [57,58]. 

The main idea of the present paper is the necessity of using representative ensembles for 
describing Bose-condensed systems. This implies that proper allowance must be made for 
all conditions which uniquely define the employed field variables. Here the consideration has 
been based on the classical Bogolubov approach [1-4] introducing two field variables, the 
condensate function r] and the operator of uncondensed particles ipi. By their definition, 
these variables are independent of each other. For a uniform system, the condensate cor- 
responds to the zero-momentum state, while this state is excluded from the description of 
uncondensed particles. This is evident from definition (86) of ■^i, where k 7^ 0. The variables 
77 and ipi are also orthogonal to each other, in agreement with definition (17), which becomes 
obvious from Eq. (86), since 



For so introduced independent orthogonal variables, it is necessary to define two normaliza- 
tion conditions and, respectively, two Lagrange multipliers controlling these normalization 
conditions. Only then there can be the assurance that the theory will be self-consistent in 
any calculations. 

Of course, if the field variables are introduced in a different way, with some other condi- 
tions, this would require to define another ensemble, with the appropriate Lagrange multipli- 
ers, whose number could also be different. For example, Hugenholtz and Pines [9], as well as 
later Gavoret and Nozieres [60] , when deriving the Hugenholtz-Pines relation on the basis of 
thermodynamic properties, did not use the standard grand ensemble. Any attentive reader 
can immediately infer from the original works [9,60] that these authors have used a different 
ensemble. They treat a uniform equilibrium system, defining the number of condensed par- 
ticles No by minimizing the internal energy at zero temperature, which corresponds to the 
minimization of free energy at finite temperatures. After defining in that way A^o — Ao(p, T), 



operator of uncondensed particles, but not of the total number of particles, as it would be 
in the standard grand Hamiltonian. Since A^o = Ao(p, T) has been explicitly substituted 
everywhere, one needs the sole Lagrange multiplier for the normalization of uncondensed 
particles. But mathematically this is absolutely equivalent to the introduction of an addi- 
tional Lagrange multiplier /iq, as is done in the present paper, before substituting A^o(P) T), 
which is defined later from the corresponding normalization condition. These ways, as is 
absolutely clear, are equivalent, but the latter method is more convenient for more general 
cases of nonuniform or nonequilibrium systems. 

Another example is given by the approach advanced by Faddeev and Popov [61] and used 
later by Popov [12-14] . They define the field operator 




(171) 



the latter is explicitly substituted into the effective Hamiltonian H — jJ^Ni, where A^i is the 




(172) 
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in which the second part in the right-hand side, 



V.^p(r) = 47 E^'^e''^^ (173) 



includes the term 0,0 7^ with k = 0. This is contrary to the Bogohibov field operator of 
uncondensed particles (86), not containing this zero-momentum term. The Fourier transform 
Cfe of the operator 

is connected with the Fourier transform of ippp by the equation 

Ck = 5ofe \/iVo + Ofc . 
For the zero-momentum state, one has 

Co = \fN^ + ao . 

Recall that in the Bogolubov case, one would have Cq = Thus, in the Faddeev-Popov 

approach, the condensate is not completely separated from uncondensed particles, but ippp 
does contain the zero-momentum term Oq. In other words, ippp is not independent from po- 
Moreover, the Faddeev-Popov representation (172) for ip consists of two parts that are not 
orthogonal to each other, 

-7^ j ^Fp(r) rfr = ao 7^ , 



V 

which is contrary to the Bogolubov case (171). Hence, ippp is neither independent of ^/po no 
orthogonal to it, but both of them define the sole variable (172). As far as, in the Faddeev- 
Popov approach, ao 7^ 0, the latter yields the interference terms in physical operators. For 
instance, here the number-of-particle operator becomes 

N = J V'^(r)V'(r) dr = No + ^4"^ + V^o (4 + «o) ■ (174) 

k 

This operator is normalized to the total number of particles, so that 

N = <7V> = No + Y.< alak > , (175) 

k 

which requires < ao >= 0. Faddeev and Popov in their original paper [61] emphasized 
that their representation (172) is principally different from the Bogolubov shift (16) and 
discussed in detail the corresponding differences. Because in the Faddeev-Popov approach 
there is a single independent variable iJj, with the normalization condition (175), the appro- 
priate representative ensemble here is the standard grand ensemble, with the grand Hamil- 
tonian H — /iN, though the number-of-particle operator (174) here is different from the 
Bogolubov form (24). More complicated forms of physical operators and the necessity to 
comply with relation (172) at each step of any calculational procedure make the usage of 
the Faddeev-Popov approach more complicated and quite inconvenient for mean-field type 
approximations. However, needless to say that a theory with one independent variable and. 



33 



respectively, with one Lagrange multiplier, is mathematically equivalent to the theory with 
two independent variables and two Lagrange multipliers. Which representation to choose is 
rather a matter of convenience. 

It is even admissible to introduce no Lagrange multipliers and to deal with the canonical 
Gibbs ensemble. But the problem with the latter is that then one has to keep the total 
number of particles fixed not merely on the average but exactly at each step of any cal- 
culational procedure. This requires to work in a restricted Fock space, where the number 
-of-particle operator degenerates to the number N = N. To accomplish this, one may resort 
to the Girardeau-Arnowitt approach [24,25] introducing the so-called number-conserving 
field operators 

This representation, however, is valid only when the number of condensed particles is large, 
A'o ^ 1, so that it is not applicable in the vicinity of the condensation point, when Nq — 0. 
Also, such a canonical representation in approximate calculations yields, as is known [24], a 
gap in the spectrum. It was noticed by Girardeau [62] and showed by Takano [63] that to 
get a self-consistent theory in the canonical ensemble requires to use the whole Hamiltonian, 
without reducing it to approximate forms. But the theory with a whole Hamiltonian for 
interacting particles has no exact solution. So, the usage of the canonical ensemble is un- 
practical in analytic investigations, though it may be employed for numerical computations 
[64,65]. 

Concluding, the choice of an ensemble is, generally speaking, a matter of convenience. 
But in any case, the chosen ensemble must be representative, which necessitates to accurately 
take into account all conditions uniquely defining the considered system. The choice of an 
ensemble is intimately connected with the choice of the appropriate field variables, which 
should not be confused. Employing inappropriate variables, that is, using a nonrcprcscnta- 
tive ensemble may result in the appearance of inconsistences and paradoxes. For example, 
resorting to the canonical ensemble, one should use the Girardeau-Arnowitt representation 
with the number-conserving field operators [24,25] or the Carusotto-Castin-Dalibard rep- 
resentation based on stochastic fields [64,65]. If one prefers the standard grand ensemble, 
then the Faddccv- Popov representation is appropriate [61]. A nonstandard variant of the 
grand ensemble, due to Hugenholtz and Pines [9], can also be used. But when the clas- 
sical Bogolubov representation [1-4] is chosen, where the field variables of condensed and 
uncondensed particles are separated from each other, being independent and orthogonal, 
then the approach developed in the present paper is the most convenient, being completely 
self-consistent. 
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Appendix A 

Here we prove relation (50). From definition (36) for tlie grand potential with Hamilto- 
nian (30), it immediately follows that 

dn on 



On the other hand, 

on 

N^- — = No + N,. 
oil 

Comparing the above equations, we get 

on 

+ 



Differentiating again the latter equation, we have 

d^n d^n d^n d^n 

+ ^ + 2 



dfi^ djjLQ dni d/iodni 
By direct calculations, we find that 

^ = -(3A\No) , ^ = -(3A\N,) , = -p cov(iVo, N,) , 

where the covariance of two operators, A and B, is 

cov(i, B) = ^ < AB + BA > - < A >< B > . 

Summarizing these equations, and using the property of the dispersion of a composite oper- 
ator, 

A2(i + 5) = A'(i) + A\B) + 2cov(i, B) , 

we obtain 

which proves relation (50). 
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Appendix B 



The fact that the terms hnear in ipi or tpl must be absent in the Hamiltonian, in order 
to satisfy constraint (25), can be proved in two ways. 

We may, first, consider a quadratic Hamiltonian approximating the exact one. Linear 
terms in such a Hamiltonian can also be present. Quadratic Hamiltonians of this type can be 
diagonalized with the help of exact canonical transformations [45]. Rigorous mathematical 
properties of the corresponding transformations, called nonuniform, are expounded in the 
book by Berezin [42] . After the Hamiltonian is diagonalized, it is straightforward to calculate 
explicitly all averages, which show that the linear terms in the Hamiltonian induce nonzero 
< V'l > a-nd, vice versa, zero hnear terms lead to < ■0i >= 0. Then one should consider 
perturbation theory, starting with the diagonalized quadratic from, and check that zero 
linear terms yield < ipi >= in all orders of the theory. This way is rather cumbersome, 
and below another method is presented. 

Let us consider the general Hamiltonian (52). Its term, linear in and is 



Mr)[-^ + uU*{r) dr 



+ j ^{r -r') ipl{r)\r]{r')\^r]{r) +r]*{r)\r]{v')\^ijji{r) drdr'. 
The equation of motion (35) can be represented as 

I — ibAr, t) = — — I H 1 . 

The first term here gives the right-hand side of Eq. (60). The second term results in the 
form 

5ipl{r,t) 

proportional to the unit operator Ijr in the Fock space generated by the field operators 

tpi, with the nonoperator complex function 



C{r,t) = 



2m 



+ U + J ^r-r')\r){r',t)\^dr' 



r]{r,t) - X{r,t) 



The equation of motion, written above, is an operator equality defined on J-'lipi). The 
operator equality assumes that it holds true at least in the weak sense implying the equality 
of all matrix elements for the states from the space the operators are defined on. In the 
present case, the latter means the space T{ipi). The vacuum state of this space is defined in 
the standard way as the vector |0 >, for which ipi{r, t)\0 >= 0. Considering for the equation 
of motion, with respect to ipi, the matrix element over the vacuum state, and, taking into 
account constraint (25), we get 

C(r,t) = 0. 

The latter equation results in the Lagrange multiplier (54), which yields H^^^ = 0. Thus, in 
order to preserve constraint (25), the linear terms in the Hamiltonian must be zero. 
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The physical meaning of the proved theorem, that < ijji >= necessarily requires 
0, is quite evident. The terms in the Hamiltonian, linear in ipi and ipl describe the physical 
processes of annihilation and creation of single particles. If such processes were permitted, 
then the function C(r, t), introduced above, is not zero. Then the evolution equation for the 
average < ^) > contains the term 



generating the nonzero value oi < ipi > . 

One should not confuse the considered situation with the often used method of specially 
adding to the Hamiltonian the terms, linear in ipi and tpl, in order to break the gauge 
symmetry. Then, of course, the average < -01 > is not zero. But at the end of calculations, 
one always set the additional linear terms to become zero, hence restoring the property 
< ipi >= 0. This procedure is what is called the Bogolubov method of infinitesimal sources 
[3,4]. These sources are usually lifted after the thermodynamic limit, but can also be made 
infinitesimally small in the process of taking the thermodynamic limit, provided this is done 
in the appropriate way [5,66]. In the approach, followed in the present paper, we do not 
need to introduce infinitesimal sources, since the gauge symmetry has already been broken 
by the Bogolubov shift (16). 
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Appendix C 



The critical temperature Tc for the Bose system with contact interactions has been ob- 
tained, in the HFB approximation, from expansions for the condensate fraction (149) and 
superfiuid fraction (150). This temperature was found to coincide with the critical tem- 
perature of the ideal Bose gas. In order to emphasize the correctness of this result, let us 
consider, first, the more general case of an arbitrary interaction potential $(r), provided 
it possesses the standard property of symmetry, such that $(— r) = $(1"), and diminishes 
sufficiently fast with increasing |r|. 

At the critical temperature when po = and pi = p, Eq. (112) reduces to 

fc^ 1 

The Fourier transform 

= j $(r)e-'(''+P) '- dv 

can be simplfied remembering that, by assumption, the interaction potential diminishes fast 
with increasing r = |r|. Then in the above integral, one can expand e~^^'^ in powers of k • r 
up to the second order. As a result, we get 

, „ ^ ( 1 - ^ khl 1 $ 



where the notation for the effective interaction radius ro is introduced, defined by the equa- 
tion 

2 _ / r^$(r)(ir 
^° " /$(r)dr ■ 

In this way, we find 

The critical temperature is given by the equation 

dk 



in which 

From here we obtain 
where the effective mass is 



(27r)3 ' 



-'-c — 



2n 



am. 



m 



2/3 



m 



1 _ r^,$,^i!^ 

J- 3 J 'ifc^A:(2^)3 



Thus, for nonlocal interactions, with an interaction radius rg, the effective mass increases, 
so that the critical temperature diminishes, as compared to the critical temperature of the 
ideal Bose gas. 
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However, for the contact interaction potential (83), we have tq = 0, hence m* = m, 
and Tc, in the frame of the HFB mean-field approximation, coincides with the ideal gas 
condensation temperature. 

This conclusion is in agreement with other studies of the critical temperature for inter- 
acting Bose gas. There exists quite a number of such investigations, as reviewed in Refs. 
[31,49]. The most accurate of these calculations are those employing Monte Carlo simula- 
tions [67-71] and those based on the optimized perturbation theory [72-74], as has been 
done in Refs. [75-80]. These investigations show that the first correction to the critical 
temperature comes from effects beyond the mean-field approximation. 
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